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Abstract 

We study thermodynamics of strongly coupled lattice QCD with two colors of massless staggered 
fermions as a function of the baryon chemical potential in 3 + 1 dimensions using a new cluster 
algorithm. We find evidence that the model undergoes a weak first order phase transition at /i = 
which becomes second order at a finite /i. Symmetry considerations suggest that the universality 
class of these phase transitions should be governed by a 0{N) x 0(2) field theory with collinear 
order, with = 3 at /U = and = 2 at /U 7^ 0. The universality class of the second order 
phase transition at ^ 7^ appears to be governed by the decoupled XY fixed point present in the 
0(2) X 0(2) field theory. Finally we show that the quantum (T = 0) phase transition as a function 
of /i is a second order mean field transition. 
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I. INTRODUCTION 



Understanding the phase diagram of Quantum Chromodynamics (QCD), as a function of 
temperature (T) and baryon chemical potential (/i) is of great interest in the phenomenology 
of strongly interacting dense matter l|. There are many excellent reviews on the subject 
and some recent ones can be found in The only first principles approach to 

the subject is based on the lattice formulation of QCD in which one computes quantities 
using the Monte Carlo method. Unfortunately, due to a variety of algorithmic difficulties 
this has been difficult to accomplish. At intermediate and large chemical potentials and 
small temperatures the numerical methods suffer from a severe sign problem. Thus, the 
most reliable lattice calculations can only be found at small u where reasonable algorithms 

riQ 

are available p, |7[. However, even these calculations become difficult especially for large 
lattices and realistic quark masses. Thus, it is fair to say that quantitatively many features 
of the T — fi phase diagram of QCD still remains unclear. A recent review of the status of 
lattice calculations at finite temperature and density can be found in 0, 1^ . 

Given the difficulties of studying the phase diagram of QCD, it is interesting to consider 
QCD-like theories which do not suffer from sign problems at /i 7^ jl^. The sign problem is 
avoided due to special properties of the fermion action which makes the fermion determinant 
real and positive. As a price, baryons turn out to be bosons. In spite of this difference QCD- 
like theories provide interesting toy models for QCD. In certain cases they are also interesting 
in their own right. A famous exam 
the years both theoretically [ll|, Il2 



e is two-color QCD and has been extensively studied, ove r 
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14j | and numerically 



20 



21 



Two-color lattice QCD with staggered fermions (2CLQCD) is especially interesting due to 
an enhanced U{2) symmetry at zero quark mass and baryon chemical potential. As we 
will see later, the long distance physics of this theory in the T — fi plane is described by 
an 0{N) X 0(2) (iV = 3 when = and = 2 when a4 0) field theory. Such field 
theories arise naturally in many condensed matter systems j23| including the studies of spin 
and charge ordering in cuprate superconductors j2J| and superfiuid transitions in ^He (25l |. 
More concretely, the normal-to-planar superfiuid phase transition in ^He is governed by an 
0(3) X 0(2) field theory 26|, which is similar to the one that arises in 2CLQCD at /x = 0. 
Although some progress has been made in uncovering important qualitative features of the 
phase diagram of 2CLQCD, many quantitative questions remain: (1) What is the order of 
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the finite temperature chiral transition at zero and non-zero chemical potentials? (2) Can the 
low energy physics at small T and fi be understood quantitatively by an appropriate chiral 



perturbation theory 



21|? (3) What is the order of the phase transition that occurs when 



the lattice gets saturated with baryons at T = 0? The reason for the lack of quantitative 
progress can be traced to the fact that the difficulty of calculations are similar to those 
in QCD at zero chemical potential. Hence, all previous studies have been limited to small 
lattice sizes and relatively large quark masses. Further, diquark condensation occurs in 
2CLQCD and is difficult to study in the conventional approach. One usually has to add a 
diquark source term and then extrapolate it to zero. 

Fortunately, new Monte Carlo algorithms have emerged for lattice gauge theories in the 



limit of infinite gauge coupling 27[ where many of the conventional algorithmic problems 
disappear. Although this strong coupling limit has the worst lattice artifacts, the qualitative 
physics is expected to be preserved. Since most of the current work is being done at couplings 
which can be considered rather strong, the phase diagram at these couplings may not be 
altered significantly as compared to the infinite coupling theory. On the other hand, thanks 
to the new algorithms, one can study the chiral limit on large lattices with ease at infinite 
coupling. Thus, studying strong coupling 2CLQCD should be a useful step in our general 
understanding of the subject. The strong coupling limit attracted the attention of physicists 



in the eighties when a variety of qualitative results were obtained using mean field theory and 
numerical work [2^ |2^ 120, iMl Interestingly, even today many 



qualitative questions continue to addressed in this limit |40, |41|, . The strong coupling 
imit of 2CLQCD was originally considered in and has been recently reviewed in 

4^. However, many interesting questions, including the ones raised above, have remained 
unanswered even in this simplified limit. 

In this article we extend the directed path algorithm invented in to study strong 
coupling 2CLQCD in the chiral limit and attempt to answer many questions including those 
raised above. Our article is organized as follows. In section |H] we discuss the model and the 
expected physics in detail. In section IIIII we explain the new algorithm which is followed 
by a section in which we discuss the observables and how we measure them. Section |3 
contains our results which is followed by a section where we present a summary of our work 
along with conclusions. A preliminary version of this work appeared in a recent conference 



proceedings 



44|. 
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II. THE MODEL 



The action of 2CLQCD we study is given by 

S = -y^Jar]a{x) 

The Grassmann valued quark fields and xi^)^ associated to the 3 + 1 dimensional 
lattice site x with coordinates {xt, xi, X2, x^), represent row and column vectors with 2 color 
components. The components will be denoted as Xa ^^"^ Xa^o, = 1,2. The gauge fields 
Ua{x) are elements of SU{2) group and live on the links between x and x + a where a = 
t, 1, 2, 3. The factor Tq, = 1 for a = 1, 2, 3 and = l/af. At weak couplings at acts as the 
temporal lattice spacing (assuming spatial lattice spacing is 1). However there is no reason 
to expect this interpretation to hold at strong couplings. Thus, we think of at as merely an 
asymmetry factor between spatial and temporal directions. It allows us to study the effects 
of temperature on asymmetric lattices and was already used for this purpose in jQ. In 
the dimer-baryon loop representation which we will construct later, the parameter l/(at)^ 
is more natural. By choosing an Lt x L'^ lattice (periodic in all directions) we can study 
thermodynamics in the L ^ oo limit at a fixed Lt by defining T = l/(at)^ as the parameter 
that represents the temperature. Zero temperature studies involve the limit ^ oo with 
fixed T. The parameter fi represents the baryon chemical potential. The absence of the 
gauge action shows that we are in the strong gauge coupling limit. 



of^t5t,^Y(^^)Ua,{x)x{x + d) - e-'"''^''''x{x + d)Ul{x)x{x) 



A. Internal Symmetries 



A detailed discussion of the symmetries of 2CLQCD can be found in |18l. I43| . For com- 
pleteness we review them briefiy here. We first rewrite eq. as 



E 



S=- > ria[x) 

X even,a=l,2,3 



Xe{x)Ua{x)Xo{x + Ct) — Xe{x)U}^{x — a)Xo{x — Ct) 



VAx) 
at 



X,{x)e'''^'''Ut{x)Xoix + i) - Xe(x)e""*'^"-^[//(a; - i)Xo{x - i) 



(2) 
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where and Xo are given by 

^=(Xe,-xrr2), Xo^i ^\ I (3) 

In our notation a arc Pauli matrices that mix x ^iid present in X„ and while r are 
Pauh matrices that act on the color space. Thus, T2UT2 = U* since U is an element of 
SU{2). 

Clearly, when /i — our model has a U{2) global symmetry: 

Xo VXo, Xe XeV\ V = exp(i(3 • a + i0) e U{2). (4) 
This symmetry is reduced to Ub{^) x Uc{^) in the presence of a chemical potential: 

C/b(1) : ^ exp(i(j30)Xo, ^ exp(-i(730) 

— — \p) 
Uc{l) : Xo exp{i(p)Xo, X^ X^ exp(-i0). 

Here Ub{1) is the baryon number symmetry — > e''^x(x), x(^) ^ xl^)^"*''^ and 

is the chiral symmetry of staggered fermions x(x) e^'^^^^^xi^): x(^) ~^ x{^)^^'^^^^^ where 

B. Properties of the Ground State 

When /i = one expects the chiral condensate, which is not invariant under the U{2) 
symmetry, to get a non-zero vacuum expectation value. Note that 

^ . ^ / T^o^2[(T^ (g) T2]Xo X E odd 

^i{x) ^ \ Z = 1,2,3 (6) 

I |Xe[crj (g) T2\a2X^ X e even 

transform as components of a three vector under the SU{2) subgroup of the U{2) symmetry 
group. It is easy to check that 

^1 = «[XiX2 + X1X2], ^2 = [-X1X2 + X1X2], ^3 = XX = [XiXi + X2X2], (7) 

Thus, the chiral condensate is the third component of a complex three vector. In addition 
all the components carry the same non-zero Uc{i) chiral charge. 

The above discussion makes it clear that the chiral condensate being non-zero is just a 
matter of choice. More generally, the ground state of the theory is such that = (n)ie*^. 



where n is some constant unit three vector. With this choice, the ground state still remains 
invariant under a U{1) subgroup given by = exp{i6 h ■ a) which implies that the U{2) 
symmetry is broken to a U{1) subgroup (note that the U{1) subgroup must be a part of 
the SU{2) subgroup of U{2)). When one says that the chiral condensate is non-zero, one 
implicitly chooses n along the third direction. This then implies ($1) = ($2) = 0. However 
when we study the effects of the chemical potential, it is natural to pick the ground state 
such that ($2) 7^ and ($1) = ($3) = which implies that the diquark condensate, 
{X1X2) = {X2X1) while the chiral condensate vanishes. Note that even though (xx) = 
the theory still breaks the Uc{l) symmetry since the diquark condensates carry a chiral 
charge. 

When /i 7^ the U{2) global symmetry is explicitly broken to Ub{1) x Uc{l)- At small n 
both the U{1) symmetries are expected to break spontaneously since the diquark condensate 
continues to be non-zero. As fi increases the density of baryons increases and at a critical 
value /ic the lattice becomes saturated with baryons which means that for fi > fic the diquark 
condensate vanishes. If this phase transition is second order, the low energy physics close 
to fic will be governed by a non-relativistic field theory. Renormalization group arguments 
indicate that this field theory is governed by mean field theory in > 2 (d represents spatial 
dimensions) .45] . 



C. Finite temperature phase transition 

At high temperatures all the symmetries are expected to be restored. This implies one 
must have a finite temperature phase transition for /i < fic- The order of this transition both 
at yU = and yU 7^ in not known. Since the order parameter at /i = is a complex 3-vector, 
its fluctuations are governed by the Landau- Ginzburg (LG) Hamiltonian of the 3-component 
complex field ipi^)- The U{2) symmetry is manifest as a 0(3) x 0(2) symmetry. Theories 
with A^-component complex scalar fields with 0{N) x 0(2) symmetry are interesting in 
condensed matter physics in describing a variety of materials j23|, and are described by the 
action (or classical Hamiltonian), 

S= [ d'^x \ d^ij* ■ d^,ip + r^* -1^ + ui^xp* ■ ipf + v\iIj ■ (8) 
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When r < 0, depending on the sign of two classes of ground states are allowed (note 
u > is necessary for stability). When f > the ground state has a spiral or helical order, 
while when f < the ground state has collinear or sinusoidal order. Since we found that 
ipi = (<l>j) = njC*^, we discover that close to the finite temperature phase transition the long 
distance physics of 2CLQCD with massless staggered fermions at /i = is described by 
the above complex field theory with = 3 and collinear order {v < 0). This field theory 
is of interest in the study of the normal-to-planar superfluid transition in ^He js^. The 
question of whether the = 3 theory with collinear order can be second order still remains 
unresolved. While the e-expansion predicts a fluctuation driven first order transition Q| 
recent results claim that a second order fixed point indeed exists 3|. As we will see later, 
in our work we find a weak first order transition. 

It is easy to argue that in the presence of the baryon chemical potential the finite temper- 
ature phase transition must be governed by a Landau- Ginsburg Hamiltonian similar to the 
one above but with N = 2. Note that the symmetry in the microscopic theory is reduced to 
f/(l) X U{1) which is manifest in the LG theory as an 0(2) x 0(2) symmetry. Further, near 
a finite temperature phase transition the presence of a chemical potential does not usually 
break charge conjugation symmetries in the low energy effective theory. For = 2 it is 
possible to define new fields 

/ 




1 
72 



1 z 
i 1 




(9) 



such that the LG Hamiltonian becomes 

2 



s 



+ 2{u + 2v)\ipi\'^\ip2\ 



(10) 



t is then obvious that when v < there is always a decoupled XY fixed point at {u + 2v) = 
4^. Using the knowledge of the critical exponents in the XY model it can be established 



that this decoupled fixed point is stable 47|, |48|. However, the fiow to this fixed point is 



rather slow so that corrections to the XY scaling can be substantial until one is very close 

n 

to the critical point ^49]. At a non-zero value of fi we indeed find a second order transition. 
A naive analysis indicates that the critical exponents are different from the XY exponents, 
however when the expected corrections to scaling are included the XY exponents can be 
used to fit our data. 
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III. DIMER-BARYONLOOP MODEL 



A. Dimer-Baryonloop Configurations 

One of the computational advantages of the strong couphng hmit is that in this hmit it 
is possible to rewrite the partition function, 

Z = J[DU][dxdx]exp{-S), (11) 

as a sum over configuration containing gauge invariant objects, namely monomers, dimers, 
and baryonloops j32, |35|, |37 



38| . Monomers are absent in the massless theory. In the case of 
2CLQCD a lattice configuration K of dimers and baryonloops is constructed as follows: (1) 
Every link of the lattice connecting the site x with the neighboring site x + a contains either 
a dimer ka{x) = 0,1,2 or a directed baryon-bond ba{x) = —1,0,1. When ka{x) = 0, it 
means that the link does not contain a dimer, while ka = 1(2) implies that the link contains 
a single (double) dimer. Similarly ba{x) = means the link does not contain a baryon-bond, 
while ba{x) = 1 means the baryon-bond is directed from x to x + a and ba{x) = — 1 means 
it is directed from x + a to x. We will also allow a to be negative. Thus, if a was positive, 
k^a{x) and 6-q(x) will represent dimers and baryon-bonds connecting x with x — a. (2) If a 
site is connected to baryon-bonds then it must have exactly one incoming baryon-bond and 
one outgoing baryon-bond. Further it cannot be connected to dimers. Thus baryon-bonds 
always form self-avoiding closed baryonloops. (3) Every lattice site x that does not contain 
a baryon-bond must satisfy the constraint 

J2K{x)=2 (12) 

a 

where the direction a in the sum takes negative values also. This implies that sites connected 
by single dimers also form a loop which we call a dimer-loop. An example of a configuration 
K is given in Figure H 
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FIG. 1: An example of a dimer-baryonloop configuration. 



B. Updating Algorithms 

Given the set {K} of dimer-baryonloop configurations described above, eq. (QJ) can be 
rewritten as 




h{x) + Mx)\)\og{T) + ^^^\y (13) 



Since the partition function is written as a sum over positive definite terms, a Monte-Carlo 
algorithm can in principle be designed to the study this system. Howe ver, the algorithm 



needs to preserve many constrains. A method to do this was developed in 



in the context 



of quantum spin models and later extended to dimer models in [27[ . We will now discuss an 
extension of these ideas to update the configurations K. In particular we consider three types 
of updates: a dimer-baryon loop flip update, a dimer update and a baryon update. Below, we 
discuss each of these updates in detail. Remember that we assume ka{x) = /c_Q,(x + d). The 
dimer update and baryon updates have been described such that this redundant information 
also gets updated automatically. 
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1. Dimer-baryon loop flip update 

This update is based on the observation that every baryon loop's orientation can be 
flipped without violating any constraints. Further a baryonloop can be converted into a 
dimer-loop and vice-versa. Thus, every loop can be in one of three states: a dimer-loop or 
a baryonloop with two different orientations. Let Y be the subset of lattice sites that are 
connected to either a dimer-loop or a baryonloop. Vy will be the number of these sites. We 
pick a site from Y at random and change the state of the loop C connected to that site to one 
of the three allowed states. The change can be accomplished using a heat-bath (or similar) 
update if we assign the following weights to the three states: a dimer-loop carries a weight 
1 while the baryonloop (in either orientation) carries the weight exp[2 nWt{C)], where 

Wt{C) = J2bt{x) (14) 

Note that Wt{C) changes sign if its orientation is flipped. 



2. The Dimer Update 

Let D be the set of sites connected by ka{x) = 1, 2 and Vd be the number of such sites. 
The dimer update changes the configuration on a subset of D. The update is as follows: 

1. A lattice site x & D is selected randomly 

2. If the site x lies on a dimer-loop, then there will be two different directions a such 
that ka{x) = 1. One of these two directions is picked at random. Else there will be 
one direction a such that ka{x) = 2. In that case this direction is picked. 

3. If the update just started and x is the first site, a virtual monomer is created at x. If 
not one is added to ka' (x) where a' is the direction from which x was reached. One is 
subtracted from ka{x) and the update moves to the neighboring site y = x + a. We 



27[ . See next 



will call X as an "active site" and ?/ as a "passive site" in the notation of 
subsection for more details. 

4. With all the neighboring sites y + a' that belong to D a non-zero weight Wa' is 
associated. If a' is a temporal direction then Wa' = T, otherwise Wa' = 1. If the 
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neighboring site does not belong to D then the weight is zero. For future reference 
we define the total dimer weight on the site y as Woiy) = J2a' Now based on 
the weights W^' a heat-bath (or similar) procedure is used to pick a new direction a'. 
One is subtracted from ka{y) and one is added to ka'{y). The update then moves to 
the neighboring site x = y + a'. 

5. If the site x is not the site which was picked in step 1, the update moves to step 2. 
Otherwise the site x will contain one virtual monomer and one direction a such that 
ka{x) = 1. With probability half the direction a is picked and the update moves 
to step 3 and with the remaining probability half the virtual monomer on the site is 
removed and the update ends. 



3. The Baryon Update 

The third update is just a minor modification of the dimer update. Let B be the set of 
sites connected to ka{x) = 2 or containing a baryonloop and Vb the number of such sites. 
The baryon update changes the configuration on a subset of B. The update is as follows: 

1. A lattice site x in i? is selected randomly. 

2. If the site x lies on a baryonloop, then there will be one direction a such that &q,(x) = 
— 1. This direction is picked. On the other hand if x is not on a baryonloop then there 
will be one direction a such that ka{x) = 2. In that case this direction is picked. 

3. If the update just started and x is the first site, a virtual "diquark" is created at x. 
Otherwise, one is subtracted from ba'{x), where a' is the direction from which x was 
reached. If ba'{x) = after the subtraction then ka' is set to 2. One is added to ba{x) 
and if ka{x) = 2 then it is set to zero. The update moves to the neighboring site 
y = X + a. We will call x as an "active site" and y as a "passive site" in the notation 



of 



. See next subsection for more details. 



4. With all the neighboring sites y + a' that belong to i? a non-zero weight Wa' is 
associated. If a' is the positive temporal direction then Wa' = Texp(2/iaf), if it is 
along the negative temporal direction then Wa' = Texp(— 2/xaj), otherwise Wa' = 1. If 
the neighboring site does not belong to B then its weight is zero. For future reference 
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we define the total baryon weight on the site y as Wsiy) = Yla' ^a'- Now based on 
the weights Wa' an over-relaxation procedure is used to pick a new direction a'. One 
is subtracted from ba{x) and if ka{x) = 2 then ka{x) is set to zero. One is added to 
ba'{x) and if ba'{x) = after the addition then ka'{x) is set to two. The update then 
moves to the site x = y + a'. 

5. If the site x is not the site which was picked in step 1, the update moves to step 2. 
Otherwise the site x will contain one virtual diquark and it would have been reached 
from the direction a'. One is subtracted from ba'{x) and if ba'{x) = after the 
subtraction then ka' is set to 2. The virtual diquark on the site is removed and the 
update ends. 

C. Active versus Passive Sites 

In the definition of the dimer and baryon updates we have defined active and passive sites. 
The passive sites play an important role during the measurement of observables. Hence we 
clarify these two class of sites further. Both the dimer and baryon updates are directed loop 
updates. They start on a site x, which is called an active site. Then they go through a 
series of sites which are referred to as passive and active alternately. Thus the second site is 
a passive site, the third is an active site and so on. If the first site is such that e{x) = 1 then 
all passive sites y, visited during the update, have e{y) = —1 and vice-versa. The weights 
Woiy) and WB{y) for passive sites encountered during the updates will play an special role 
in the measurement of correlation functions as discussed below. 



D. Detailed Balance and Ergodicity 

Each of the three updates satisfy detailed balance. The proof of detailed balance for 
the the dimer-baryon loop flip update is straight forward and conventional. On the other 
hand the proofs for the dimer and baryon updates need some understanding of directed 
loop algorithms. Once this is clear, the proof is essentially straight forward. We refer the 
reader to 27, and do not prove detailed balance of these two algorithms in this article. 
The combination of the three updates makes the algorithm ergodic. To see this we note 
that there is always possible to flip all the baryonloops into dimer-loops. Once this is done 



12 



one can use the proof given in [27^ to show that dimer updates are ergodic in the space of 
configurations that purely consist of dimers. 

IV. OBSERVABLES 

A variety of observables can be measured with our algorithm. We will focus on the 
following: 

(a) The chiral two point function is given by 

Gc{z,z') = [x{z)x{z)x{z)x{z)) (15) 
and the chiral susceptibility xc is 



^Y.Gc{z.z') (16) 



z' 

where V = L^Lf. Both these observables can be measured easily during the dimer 
update. It is possible to show that [27 1 

where x is the first site of the dimer update and the sum is over all the passive sites y 
encountered during the update. Vd and Woiy) were defined in the previous section. 

(b) The diquark two point function is given by 

Gsiz, z') = (xi{z)X2iz)x2{^')xi{z')^ (18) 
and the diquark susceptibility is given by 

z' 

These can be computed during the baryon update by using 

(20) 



Gb{z,z') = /j2 
\ y 



VBSx,zSy,z' 



where x is the first site picked during the update and the sum is over all the passive 
sites y encountered during the baryon update. Vb and WB^y) were defined in the 
previous section. 
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(c) Baryon density ub is defined as 



1 dlnZ , , 



and in the dimer-baryon loop flip update it can be measured by the formula 

_IVy Wm sinli(2//U;(C)) \ 

\ y ^(C) l + 2cosh(2AtW^t(C))/ ^ ^ 

where C is the dimer-baryon loop picked, Vy jV is the fraction of the lattice sites that 
contain dimer-baryon loops, 5'(C) is the size of the loop and (defined earlier) is 

the temporal winding of the baryon loop. 

(d) The helicity modulus associated with the U{\) chiral symmetry Yc 

1 ^11^ \ 

Yc = 



V E ((T.M^)] ) (23) 

* M=l,2,3 \ \ a; / / 



where 

)) (24) 

and Vs = L\ 

(e) The helicity modulus associated with the U{1) baryon number symmetry Yb 

Yb 



V E (fE^.(-)) ) (25) 

* //=1,2,3 \\x J I 



3K 

where 

B,{x) = [h,{x)] (26) 

Both Yc and Yb are obscrvables that can be calculated easily for each configuration K 
and averaged. Note also that our definitions of Yc and Yb are natural at finite temperatures. 

V. RESULTS 

A. Zero Chemical Potential 



In a finite volume there is no spontaneous symmetry breaking. However, the effects of 
symmetry breaking can still be studied by examining the large volume limit of xc ^-nd Xb] 
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if the symmetry is broken in the infinite volume limit, we expect these susceptibilities to 
grow with the volume of the system. From symmetry considerations at = it is possible 
to show that Gc{z,z') = 2Gb{z^z') which implies xc = '^Xb- Symmetry breaking can also 
be observed through the helicity modulus Yq and Yg; both must reach a non-zero constant 
if the symmetry breaking pattern is as expected. All these can be understood quantitatively 
using the low energy effective action 

^eff = J d'^x {^{d,S) ■ (d^S) + ^(d^u) ■ id,u)} (27) 

where S{x) and u{x) are unit three and two vector fields respectively. Finite size scaling 
formula for various quantities can be obtained following the discussion in ^|. We note that 
this approach to low energy physics is equivalent to others found in the literature 10|, |ll|, |21 1 . 

At a fixed value of Lf the parameter T can be increased to induce a phase transition 
between a low temperature phase with spontaneous symmetry breaking and a high temper- 
ature symmetric phase. In order to study this phase transition we have performed extensive 
calculations at a fixed Lj = 4 for different spatial lattice sizes L varying from 16 to 256 and 
for many different values of T. The low energy effective theory introduced in eq. ()27j) with 
d = 3 can then be used to predict the signatures of the broken phase. We have studied 
two such signatures: (1) Yc and Yr go to non-zero constants at large L. Extending the 
calculations of one can show 5^ 

r.2 b c 2B^ b' c' 

(2) The finite size scaling of the chiral susceptibility is given by 52| 

Xc = ^{L' + P^{^, + j,)L^}+aL (29) 

where jSi = 0.226 is the shape coefficient for cubic boxes |51] and S/v^ = (xx)- Figure El 
gives our results at T = 2.918 which is a value of T in the broken phase. The graph shows 
that the above expectations are satisfied well. 

Figure El shows the dependence of xc as a function of L for different temperatures. Using 
the data for L < 256 we find that xc increases as for T < 2.9275, but saturates for 
T > 2.9285 as L becomes large. Thus, we think is between these two temperatures. For 
a second order transition, close to Tc, one expects 

Xc^"-' = [go + 9i{T/n - 1)1'/" + ...], (30a) 
YcL = [fo + MT/n - l)L'^' + ....], (30b) 
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FIG. 2: The inset shows the data for Yc and Yg as a function of L. The hnes shown are fits to 
eq. (pS)) which yields = 0.0965(5) and i?^ = 0.0839(6). The main figure shows the pfot of xc 
versus L for T = 2.918. The soHd hne is a fit to eq. H29|) with and fixed to the values quoted 
above. The fit yields S = 0.3364(7) and a = 2.6(1) with a x^/DOF = 0.5. 



which was recently observed in other strong coupling theories |53L |5J, |55[ . Our data does 
not fit well to this form. Clearly, the data for xc shows more structure than can be captured 
by the above relations. We have also verified that YcL does not seem to scale as a constant 
for large L anywhere in the range T = 2.928 ± 0.002. Hence we think that the transition 
is not second order. Interestingly, we are able to fit the non-monotonic behavior in xc at 
T = 2.9285 and 2.9292 using the relation 

_ a + 6L3exp(-A^L3) 

l + cexp(-A^L3) ^ 

as long as we use data for L > 48. This form is natural in the presence of two phases 
(one broken and one symmetric) whose free energy densities differ by AjF. This leads us to 
conclude that the phase transition is indeed first order. 

The existence of a first order transition implies that the correlation lengths at the critical 
point do not diverge. In that case how big are these lengths at the transition? In the high 
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FIG. 3: This figure shows the plot of xc versus L for different values of T across the phase 
transition. The solid lines for T < 2.9275 are fits to eq. (|29() while those for T = 2.9285 and 
T = 2.9292 are fits to eq. (EH). We find d = 1650(200), 1150(50), b = 0.030(5), 0.027(4), c = 
2.5(5), 2.0(4) and AF = 2 x 10"'^,4 x 10"*^ for the two temperatures. Ah the fits have x^/DOF 
less than 1. 

temperature phase one can compute screening masses M^-, Mb from the exponential decay 
of Gc{z,z') and Gb{z,z') respectively, for large spatial separations between z and z'. At 
fi = we expect = Mb- In the broken phase F"^ and have dimensions of mass and 
are the relevant physical scales in the problem. In figured we show the behavior of Mb, 
F"^ and which have been obtained after extrapolations to infinite volumes. As can be 
seen, the correlation lengths close to the transition are extremely large, about 40 — 50 lattice 
units, indicating that the transition is a weak order transition. If the transition was second 
order, we would have expected = ao(^c — Ty, B^ = h^iTc — T)" and Mb = Co(T — T^Y . 
Interestingly, these relations describe the data reasonably well, but with different and 
V in the two phases. In the low temperature phase a combined fit of F"^ and B^ gives 
ao = 0.562(3), 6o = 0.488(3), = 2.92856(4) and u = 0.387(2) with a x^/DOF = 1.5. On 
the other hand the fit of Mb gives cq = 0.64(5), = 2.9266(6) and u = 0.50(3) with a 
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n l_i I I I I I I I I I I I I I I J J I I I I I I I I I I I I I I I 

2.9 2.91 2.92 2.93 2.94 2.95 2.96 

T 

FIG. 4: This figure shows the plot of Mb, F"^ and B'^ as a function of T. The sohd hues are fits 
to the data as discussed in the text. 

/ DOF = 2.2. A combined fit of all the data on both sides of the transition with a single 
Tc and u does not fit well confirming our claim that the transition is first order. 

B. Non-zero Chemical Potential 

Having established a first order transition at zero chemical potential we next focus on the 
finite temperature transition at /i = 0.3 with Lt = 4. The chemical potential reduces the 
symmetry of the theory to Ub{1) x Uc{1)- At low temperatures both the U{1) symmetries 
are broken leading to two Goldstone bosons. The two correlators Gc{z, z') and Gb{z, z') are 
no longer related: Gc{z,z') decays exponentially while Gb{z,z') remains non-zero for large 
\z — z'\. This means xc saturates for large L while xb grows with the volume and shows 
the presence of a diquark condensation and baryon superfluidity. Yc and Yb again go to 
non-zero constants at large L. In order to determine the finite size scaling formula, we use 
the same effective theory as given in eq. (|77j) except that now both S{x) and u{x) are unit 
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two vectors. The large L limits of Yc and are now given by 

'^c = F' + j + ^ + .., Ys = B' + j + ^,+.... (32) 
The effective field theory also predicts that xb is given by 

XB = ^{L' + f3,{j-^ + ^,)L'] + aL. (33) 

where A/a/I^^ = {X1X2) = {X2Xi)- 

Interestingly, we find that the values of Yb and Yc come together as we increase the 
chemical potential, although for small temperatures and small chemical potentials we can 
still distinguish between them. On the other hand close to the finite temperature phase 
transition they become indistinguishable. We note that the action in eq. (fTUj) predicts 
F"^ = close to the phase transition. The operator which splits them is irrelevant and goes 
to zero at the transition. However, this explanation does not explain why Yc and Yb come 
close to each other as a function of /i. This behavior should be examined using effective 
field theory and may emerge naturally, but we have not attempted it so far. Figure El gives 
our results at T = 2.85, a value of T in the broken phase. The inset shows the behavior of 
iYc looks identical within errors). Fitting the data, we find = = 0.0378(11) with 
a x^/DOF of around 0.5. Fixing F"^ and B^ to these values, our data for xb fits "well to 
eq. dSni) as long as we use L > 24. We get A = 0.117(1) and a = 2.4(3) with x^/DOF = 1. 
The fit is shown as a solid line in the figure. 

Figure IHl also shows at T = 2.86 and 2.87 as a function of L. Unlike the /i = case, 
Xb behaves monotonically suggesting that the transition could be second order. In order to 
check this we can verify if eq. ()30b|) is satisfied close to Tc. A fit of our data to this relation 
gives Tc = 2.85946(7), u = 0.60(2) /o = 1.128(3) and /i = -0.161(3) with a x^/DOF = 1. 
In figure El (plot on the left), we show our data and the fit. In the inset of the figure we plot 

at T = 2.8594. Since this value of T is close to Tc we expect xb = doL^^"^ should describe 
the data reasonably well. Indeed a fit shows go = 0.306(8), r] = 0.042(2) with a x^/DOF of 
1.2 and is shown as a solid line in the inset. The values of T^, A^ and Mb obtained from 
the infinite volume extrapolations must also satisfy 

F^ = ao(Tc - T)^ Mb = co{T - T,)u, A^ = do{T, - Tf". (34) 

Figure El (plot on the right) shows a combined fit of these three quantities as a function of 
T using the value of T^ obtained above. We find v = 0.610(6), f3 = 0.311(5), Oq = 0.70(1), 
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FIG. 5: This figure shows XB as a function of L at /i = 0.3 at three different temperatures close 
to the critical temperature. The inset shows that Yb goes to a constant as a function of L at 
T = 2.85. The solid line is a plot of eq. ^ with A = 0.117 and a = 2.4(3). Clearly, xb grows 
as at T = 2.85 but be gins to saturate at T = 2.87 indicating that there is a transition between 
these two temperatures. 

Co = 0.79(1) and do = 0.252(3) with a x^/DOF = 1.3. Note that these critical exponents 
do satisfy the hyper-scahng relation (3 = {d — 2 + vi)u /2. Thus, our data strongly supports 
the existence of a second order transition. 

Unfortunately, the above results are in contradiction with the expectation from section 
III CI where it was argued that the critical behavior at 7^ must be governed by three 
dimensional XY universality class. This implies that we should have obtained v = 0.6715, 
/3 = 0.3485 and r/ = 0.0380 Q and not the exponents we found above. As discussed in 
section Hi C\ the problem is that in a model with U{1) x f/(l) symmetry and coUinear order, 
the corrections to the XY scaling, due to the irrelevant direction in the u, v plane (see section 
III C|) . are rather large. Taking into account the leading corrections to scaling one expects 



YcL = fo + f^L-^ + fi{T- Tc)L^I^ close to where uj = 0.0218 |49j. The smallness of uj 



makes the corrections large for the lattice sizes we have explored. In fact we find that our 
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FIG. 6: (Left) Plot of YcL versus T close to Tc for various values of L. The solid lines are plots 
of YcL = 1.128 - 0.161(r - 2.85946)Li/°-6°. The inset shows xb as a function of L at T = 2.8594 
and the solid line is a plot of CSOeL^-^ss. (Right) Plot of Mb, F^A^ versus T close to Tc. The 
solid lines show the plots of eq. (EH) with oq = 0.70, cq = 0.79,(io = 0.252, u = 0.610, (3 = 0.311 
and Tc = 2.85946. 



data also fits well to this corrected form if we use the XY critical exponents and the known 
value of uj in the fits. A combined four parameter fit of our data close to Tc now yields 
Tc = 2.8590(3), /o = 2.321(2), = -1.260(1), /i = -0.299(6) with a x^/DOF = 1.6. This 
fit is shown in figure [3 (the plot on the left). At T = 2.8590, xb fits well to the form goL"^'^ 
when T] is fixed to 0.0380. The fit gives = 0.302(4) with a x^/DOF of 0.57 (solid line in 
the inset of figure [Tj). When the scaling corrections are included in F"^, and Mb one gets 



F = ao(T, - TYil + a'o(T, - T^, (35a) 
Mb = co(T - TcYil + c'oiT - TcD, (35b) 
= do{Tc - Tf^il + d',iTc - TD- (35c) 

Fixing Tc = 2.859 and using the XY critical exponents and u as above, a combined fit again 
works very well and is shown in the plot on the right in figure [7| We get = 4.8(5), = 
-0.8(1), Co = 4.6(4), c[) = -0.83(8), do = 2.1(1) and d'^ = -0.89(5) with a x^/DOF = 1. 
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FIG. 7: Same as figure IHl but now showing the solid lines show 0(2) scaling including corrections. 
(Left) The solid lines represent YcL = 2.321 - 1.26L-o °2i8 - 0.299(r - 2.8590)L^/°-^'^^"^ . The inset 
shows xb as a function of L at T = 2.8590 and the solid line is a plot of 0.302L^'^^^. (Right) 
The solid lines are plots to eqs.(jSni) with ao = 4.8(5), a'o = -0.8(1), cq = 4.6(4), Cq = -0.83(8), 
do = 2.1(1), d'o = -0.89(5), u = 0.6715, /3 = 0.3485, w = 0.0218. 

C. Zero Temperature 

Next we turn to the physics at zero temperature. For this purpose we compute quantities 
with Lf = L at T = 1.0 for various values of fi and L. We now expect 



Xb 



IS? 
2 



(36) 



where A = {X1X2) = {X2X1) 7^ 0. Note that since we use the same finite size scaling form 
in four and three dimensions, A is normalized differently here as compared to the finite 
temperature case. The chiral susceptibility xcy must again saturate with L at any non-zero 
value of yU. Finally, due to our definitions (see eqs. pH|l . (j^ ) both the helicity modulus 
Yc and Yb, grow linearly with L for large L. These expectations emerge nicely in our 
calculations as can be seen in figure |H1 

As the chemical potential increases the average number of baryons in the ground state 
increases. At some critical chemical potential fic, the ground state has a baryon on every 
lattice site. Since the baryons behave as hard-core bosons, at fi^ there must be a phase 
transition to a phase where superfluidity is no longer present. We now focus on this phase 
transition. Renormalization gro up arguments show that this phase transition must be gov- 
erned by mean field theory A mean field analysis was performed recently in 4^ and 
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FIG. 8: Results at T = 1.0 and /i = 0.01. Note that xb grows with the volume while xc saturates. 
Both Yq and Yb grow linearly with L. The lines are drawn to guide the eye. 

the critical chemical potential was found to be fi^ = 0.5 cosh"^(v^) = 0.909223... The 
diquark condensate was shown to be 




In order to check these results we have computed the diquark condensate A by fitting our 
Xb data to the relation xb = ^[L'^ + + b"]. Figure IHl shows our data along with the 
mean field result j3] and the result with one- loop corrections j^^. Clearly, the one- loop 
corrections are necessary before connection with mean field theory can be established. 

The fact that this phase transition is driven due to the saturation of the lattice with 
baryons can be seen in the inset of figure where the baryon density is plotted as a function 
of /i. For fi > fic it costs an energy Eb to remove a baryon. This energy gap grows linearly 
as {fi — fic) ■ Thus for fi > fic one obtains a phase containing non-relativistic particles whose 
dispersion relation for small momenta looks like E{p) = Eb+p"^ /2Mk which leads to a spatial 
correlation length ^ = 1/ \/2EBMk- Since we expect ^ to scale as 1 / \/{fi — fic) close to fic 
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FIG. 9: The figure shows our results for A, and Eb as a function of fi. For fj, < the 
dotted line is the mean field result for A given in eq. (|37|) and the solid line contains the one loop 
corrections. For fj. > fic the dashed line is the plot of 3.549^/x — fic and the solid line is the plot of 
2(/i — fic)- We have used fic = 0.90922 here. The inset shows the baryon density ns as a function 
of fl. 

one expects the kinetic mass Mk to be a constant. Figure IHl shows the plot of Eg and 
as a function of fi. We find ^"^ = 3.549(4) ^(/x - 0.90920(3)) with a x^/DOF = 0.26 and 
Eb = 2.000(2) (/i - 0.90922(1)) with a x^/DOF = 0.35. Again /i^ is in excellent agreement 
with the mean field result which provides strong evidence that the phase transition is second 
order and belongs to the mean field universality class. From the behavior of Eb and ^ we 
find Mk = 3.12(5). 

VI. DISCUSSION AND CONCLUSIONS 

In this work we constructed an efficient cluster algorithm and studied the phase structure 
of two color lattice QCD with massless staggered fermions in the strong coupling limit. We 
found that the finite temperature phase transition at zero chemical potential is weakly first 
order, while the same transition at an intermediate value of the chemical potential was 
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FIG. 10: The phase diagram of two-color lattice QCD with staggered fermions at strong couplings. 

second order. This second order transition was found to be in the universahty class of the 
three dimensional XY model as expected from theoretical arguments. However, in order 
to show this we needed to include the large corrections to scaling expected in the theory. 
The quantum phase transition at zero temperature between a baryon superfiuid phase and 
a normal phase was also found to be second order in the mean field universality class. The 
physics in the normal phase was that of interacting non-relativistic particles. Based on 
these observations, we can attempt to draw the full diagram of two color QCD in the strong 
coupling limit. Our proposal is shown in figure ITUl 

It is important to understand how this phase diagram will change as we go to weaker 
couplings and towards the continuum limit. Clearly, when the four flavor nature of staggered 
fermions becomes important, it can begin to change significantly. However, for couplings 
that are currently explored, universality arguments suggest that the phase diagram will 
qualitatively remain the same although quantitatively the values of the critical points will 
mostlikely get affected. It is also interesting to ask why the transition at /i = is so weak. 
As discussed in section [Hi renormalization group studies based on the e-expansion indicate 
that the transition will generically be a fluctuation driven first order transition [23^. Such a 
transition could be weak as we observe. On the other hand, recent work based on high order 
perturbation theory combined with resummation techniques do not rule out the possibility 
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of a second order transition [26[ . This means we are in the vicinity of a tricritical point and 
by changing some parameter in the theory we could change the transition to second order. 
When this occurs, it is hkely that the weakly first order line in the above phase diagram 
will disappear. Such a scenario can in principle be investigated by introducing more tunable 
parameters within our model and by studying their effects on the order of the transition. 

It is amusing to note that the above phase diagram is different from the standard phase 
diagram in QCD where the baryon chemical potential induces a first order transition instead 
of weakening it. This non-standard scenario was discussed in as a possibility in QCD and 
has also been seen in the potts model simulations Finally, we note that the physics close 
to the quantum critical point is also interesting from the point of view of non-relativistic 
field theory. Since the model has a U{1) x U{1) symmetry, the low energy effective field 
theory here is richer than in a theory with a simple U{1) particle number symmetry. 
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Appendix A: Exact results on a 2 x 2 Lattice 



Below we list the exact expressions for various observables on a 2 x 2 lattice. We have 
tested the algorithm against these exact results. Tables E] and in] show the comparison of 
exact results against those obtained using the algorithm. 



1 8T3(3 + 2cosh(^)) + 32(T2 + T) +40 

'4 9T4 + 2T4(1 + cosh(^) + 6 cosh(^)) + 25 + IQT\2 + cosh(^)) 



1 (48T3 + 32T3 cosh(^) + 64T) cosh(^) + 8(1 + cosh(i|))4T2 + 80 
8 9T4 + 2T\l + cosh(^) + 6 cosh(i|)) + 25 + IQT\2 + cosh(i|)) 



1 8T4 sinh(||) + 16(2T2 + l.^T^) sinh(i|) 

'''' " 4 9T4 + 2T\l + cosh(^) + 6 cosh(^)) + 25 + IQT\2 + cosh(^)) 

1 32T2 cosh(^)) + 32T2 + 80 

" 4 9T4 + 2T4(1 + cosh(^) + 6cosh(i|)) + 25 + IQT^{2 + cosh(i|)) ^^^^ 

1 32T2 cosh(4fe)) + 64T2 + 80 

c 4 9T4 + 2T4(l + cosh(^) + 6cosh(i|)) + 25 + 16T2(2 + cosh(i|)) ^ ' 
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XC 


Xb 






exact 


algorithm 


exact 


algorithm 


exact 


algorithm 


0.0 


0.734694 


0.73463(33) 


0.367347 


0.36740(6) 


0.0 


0.0 


0.1 


0.719686 


0.71926(27) 


0.366698 


0.36674(6) 


0.074563 


0.07441(8) 


0.5 


0.412785 


0.41267(31) 


0.324066 


0.32405(7) 


0.462172 


0.46224(20) 


2.0 


0.001343 


0.00138(3) 


0.018948 


0.01895(1) 


0.997655 


0.99764(3) 



0.0 



0.1 



0.5 



2.0 



Yb 



('xa.c( 



0.367347 



0.363055 



0.254861 



0.001339 



ilfforitlim 



0.36740(26) 



0.36347(21) 



0.25487(20) 



0.00135(2) 



Yc 



exact 



0.448980 



0.442306 



0.289955 



0.001340 



algorithm 



0.44887(34) 



0.44258(32) 



0.28981(24) 



0.001337(13) 



TABLE I: Exact versus Monte Carlo results at T = 1. 





XC 


XB 


riB 




exact 


algorithm 


exact 


algorithm 


exact 


algorithm 


0.0 


0.302981 


0.30305(11) 


0.151491 


0.15153(10) 


0.0 


0.0 


0.1 


0.299581 


0.29958(11) 


0.150955 


0.15095(4) 


0.082043 


0.08218(13) 


0.5 


0.229674 


0.22987(12) 


0.137719 


0.13763(2) 


0.403837 


0.40373(17) 


2.0 


0.012903 


0.01295(5) 


0.034431 


0.03444(1) 


0.966865 


0.96678(9) 





Yb 


Yc 




exact 


algorithm 


exact 


algorithm 


0.0 


0.066076 


0.066028(55) 


0.095085 


0.095026(87) 


0.1 


0.065598 


0.065625(56) 


0.094062 


0.094109(66) 


0.5 


0.054750 


0.054770(66) 


0.072868 


0.072927(78) 


2.0 


0.004203 


0.004195(13) 


0.004284 


0.004277(13) 



TABLE II: Exact versus Monte Carlo results at T = 3. 
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